#include"system_gammow.h"
#include<cmath>
#include<algorithm>
#include"help.h"
using std::abs;
using std::min;
using std::max;

void System_Tz_Gammow::setupOrbitals()
{
  spsolver.init();
  spsolver.solve();
  //  int numOfOrbitals=spsolver.Orbitals.size();
  int numOfOrbitals=50; 
  spsolver.setupOverlap(numOfOrbitals);
  Orbitals.resize(numOfOrbitals);
  for(int i=0;i<numOfOrbitals;i++)
    Orbitals[i]=spsolver.Orbitals[i];
  setupGroups();
  FermiSurface=spsolver.FermiSurface;
  HFenergy=spsolver.Energy;
  printOrbitals();
}

/*
  dim of MatEle=1,the 1 double are:
      0  interaction
*/
const int System_Tz_Gammow::MatEledim=4;
void System_Tz_Gammow::setupTwoBodyMat()
{
  setupIndiceIn2BStates();
  int totalChannels=TwoBodyChannels.size();
  TwoBodyMat.resize(totalChannels);
  for(int channel=0;channel<totalChannels;channel++)
    {
      int num2BInChannel=TwoBodyChannels[channel].size();
      TwoBodyMat[channel].clear();
      TwoBodyMat[channel].resize( (num2BInChannel*(num2BInChannel+1))/2 );
      for(int bra=0;bra<num2BInChannel;bra++)
	for(int ket=bra;ket<num2BInChannel;ket++)
	  {
	    int ab=TwoBodyChannels[channel][bra];
	    int cd=TwoBodyChannels[channel][ket];
	    int a=TwoBodyStates[ab].a;
	    int b=TwoBodyStates[ab].b;
	    int c=TwoBodyStates[cd].a;
	    int d=TwoBodyStates[cd].b;
	    int J=TwoBodyStates[ab].J;
	    TwoBodyMatEle MatEle(MatEledim);
	    get2BmatOnNeed(a,b,c,d,J,MatEle);
	    set2Bmat(channel,bra,ket,MatEle);
	  }
    }
}

void System_Tz_Gammow::get2BmatOnNeed(int a,int b,int c,int d,int J,TwoBodyMatEle & MatEle) const
{
  MatEle.resize(MatEledim);
  MatEle.setZero();
  int jja=Orbitals[a].jj;
  int jjb=Orbitals[b].jj;
  int jjc=Orbitals[c].jj;
  int jjd=Orbitals[d].jj;
  int Tzab=(Orbitals[a].tz+Orbitals[b].tz)/2;
  int Tzcd=(Orbitals[c].tz+Orbitals[d].tz)/2;
  if(a==b && J%2) return;
  if(c==d && J%2) return;
  if(Tzab!=Tzcd) return;//isospin projection check
  if((Orbitals[a].l+Orbitals[b].l+Orbitals[c].l+Orbitals[d].l)%2) return;//parity check
  int Jmin=max(abs(jja-jjb)/2,abs(jjc-jjd)/2);
  int Jmax=min((jja+jjb)/2,(jjc+jjd)/2);
  if(Jmin>Jmax) return;//angular momenta check
  
  for(int ia=0;ia<Orbitals[a].Coeff.size();ia++)
    for(int ib=0;ib<Orbitals[b].Coeff.size();ib++)
      for(int ic=0;ic<Orbitals[c].Coeff.size();ic++)
  	for(int id=0;id<Orbitals[d].Coeff.size();id++)
  	  {
  	    System_TzType::TwoBodyMatEle Temp;
  	    int iia=pSystem_Tz->Groups[Orbitals[a].GroupIndex][ia];
  	    int iib=pSystem_Tz->Groups[Orbitals[b].GroupIndex][ib];
  	    int iic=pSystem_Tz->Groups[Orbitals[c].GroupIndex][ic];
  	    int iid=pSystem_Tz->Groups[Orbitals[d].GroupIndex][id];
  	    pSystem_Tz->get2BmatAt(iia,iib,iic,iid,J,Temp);
  	    if(Temp.empty()) continue;
  	    if(iia==iib) Temp*=sqrt(2.);
  	    if(iic==iid) Temp*=sqrt(2.);
	    MatEle+=TwoBodyMatElement<complexd>(Temp)*(Orbitals[a].Coeff[ia]*Orbitals[b].Coeff[ib]*Orbitals[c].Coeff[ic]*Orbitals[d].Coeff[id]);	    
	    //  	    MatEle[0]+=pSystem_Tz->get2B(realA,iia,iib,iic,iid,J)*(Orbitals[a].Coeff[ia]*Orbitals[b].Coeff[ib]*Orbitals[c].Coeff[ic]*Orbitals[d].Coeff[id]);
  	  }
  if(a==b) MatEle*=sqrt(2.)/2;
  if(c==d) MatEle*=sqrt(2.)/2;
}

KDataType System_Tz_Gammow::get1B(int bra,int ket) const
{
  KDataType val=0;
  //  if(Orbitals[bra].isInAGroupWith(Orbitals[ket]) )
  if(Orbitals[bra].channel==Orbitals[ket].channel)
    {
      for(int i=0;i<Channels[Orbitals[bra].channel].numPoints;i++)
	{
	  KDataType k=Channels[Orbitals[bra].channel].k[i];
	  //	  KDataType wk=Channels[Orbitals[bra].channel].wk[i];
	  val+=Orbitals[bra].wavefun(i)*Orbitals[ket].wavefun(i)*k*k;
	}
      val*=(1.-1./realA)*hbar_c*hbar_c*0.5/mc2;
    }
  return val;
}

KDataType System_Tz_Gammow::get1B_p2(int bra,int ket) const
{
  KDataType val=0;
  //  if(Orbitals[bra].isInAGroupWith(Orbitals[ket]) )
  if(Orbitals[bra].channel==Orbitals[ket].channel)
    {
      for(int i=0;i<Channels[Orbitals[bra].channel].numPoints;i++)
	{
	  KDataType k=Channels[Orbitals[bra].channel].k[i];
	  //	  KDataType wk=Channels[Orbitals[bra].channel].wk[i];
	  val+=Orbitals[bra].wavefun(i)*Orbitals[ket].wavefun(i)*k*k;
	}
      val*=hbar_c*hbar_c*0.5/mc2;
    }
  return val;
}


void System_Tz_Gammow::print1B(const std::string outfile) const
{
  std::ofstream fout(outfile.c_str());
  int numOrbitals=Orbitals.size();
  for(int i=0;i<numOrbitals;i++)
    for(int j=0;j<numOrbitals;j++)
      {
	fout<<i<<"\t"<<j<<"\t"<<get1B(i,j)<<std::endl;
      }
}

void System_Tz_Gammow::print2B(const std::string outfile) const
{
  std::ofstream fout(outfile.c_str());
  int totalChannels=TwoBodyChannels.size();
  for(int channel=0;channel<totalChannels;channel++)
    {
      int num2BInChannel=TwoBodyChannels[channel].size();
      for(int bra=0;bra<num2BInChannel;bra++)
	for(int ket=bra;ket<num2BInChannel;ket++)
	  {
	    int ab=TwoBodyChannels[channel][bra];
	    int cd=TwoBodyChannels[channel][ket];
	    int a=TwoBodyStates[ab].a;
	    int b=TwoBodyStates[ab].b;
	    int c=TwoBodyStates[cd].a;
	    int d=TwoBodyStates[cd].b;
	    int J=TwoBodyStates[ab].J;
	    TwoBodyMatEle MatEle(MatEledim);
	    get2BmatOnNeed(a,b,c,d,J,MatEle);
	    //	    get2Bmat(channel,bra,ket,MatEle);
	    fout<<a<<"\t"<<b<<"\t"<<c<<"\t"<<d<<"\t"<<J<<"\t"<<MatEle[0]<<std::endl;
	  }
    }
     
}

void System_Tz_Gammow::printOrbitals() const
{
  for(int i=0;i<Orbitals.size();i++)
    {
      std::cout<<i<<"\t"<<Orbitals[i].l<<"\t"<<Orbitals[i].jj<<"\t"<<Orbitals[i].tz<<"\t"<<Orbitals[i].e<<std::endl;
    }
}

// KDataType System_Tz_Gammow::rL(int bra,int ket,int L) const
// {
//   KDataType val=0;
//   for(int ibra=0;ibra<Orbitals[bra].Coeff.size();ibra++)
//     for(int iket=0;iket<Orbitals[ket].Coeff.size();iket++)
//       {
// 	int iibra=pSystem_Tz->Groups[Orbitals[bra].GroupIndex][ibra];
// 	int iiket=pSystem_Tz->Groups[Orbitals[ket].GroupIndex][iket];
// 	val+=pSystem_Tz->rL(iibra,iiket,L)*Orbitals[bra].Coeff[ibra]*Orbitals[ket].Coeff[iket];
//       }
//   return val;
// }
// KDataType System_Tz_Gammow::Q(int p,int h,int L,int isospin) const
// {
//   int Tzh=Orbitals[h].tz;
//   int Tzp=Orbitals[p].tz;
//   if(Tzh!=Tzp) return 0;  
//   int jjp=Orbitals[p].jj;
//   int jjh=Orbitals[h].jj;
//   int lp=Orbitals[p].l;
//   int lh=Orbitals[h].l;
//   KDataType temp,fac;
//   if(L==0)
//     {
//       temp=rL(p,h,2);
//       temp*=reducedYMatEle(lp,jjp,lh,jjh,0)* sqrt(4*Pi);
//     }
//   else
//     {
//       temp=rL(p,h,L);
//       temp*=reducedYMatEle(lp,jjp,lh,jjh,L);
//     }
//   fac=isospin?Tzh:1.0;
//   temp*=fac;
//   return temp;
// }

// KDataType System_Tz_Gammow::Q1_v(int p,int h) const
// {
//   int Tzh=Orbitals[h].tz;
//   int Tzp=Orbitals[p].tz;
//   if(Tzh!=Tzp) return 0;
//   KDataType temp=rL(p,h,1);
//   int jjp=Orbitals[p].jj;
//   int jjh=Orbitals[h].jj;
//   int lp=Orbitals[p].l;
//   int lh=Orbitals[h].l;
//   temp*=reducedYMatEle(lp,jjp,lh,jjh,1);
//   temp*=(Tzh*0.5-0.5+Z*1.0/A);
//   return temp;
// }
// KDataType System_Tz_Gammow::Q1_s(int p,int h,KDataType Rms) const
// {
//   int Tzh=Orbitals[h].tz;
//   int Tzp=Orbitals[p].tz;
//   if(Tzh!=Tzp) return 0;
//   KDataType temp=rL(p,h,3)-rL(p,h,1)*Rms*5./3;

//   int jjp=Orbitals[p].jj;
//   int jjh=Orbitals[h].jj;
//   int lp=Orbitals[p].l;
//   int lh=Orbitals[h].l;
//   temp*=reducedYMatEle(lp,jjp,lh,jjh,1);
//   return temp;
// }
